!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This file contains the main program that calls the 
! other modules and subroutines as necessary, and is used 
! in estimation.
! 
! Refer to solution.f90 for the algorithms that are used
! to solve for the stationary general equilibrium of the
! economy, to simulate a panel of firms, to construct
! the model counterparts of targeted moments and other
! quantities of interest (e.g., the levels of output,
! consumption, etc.), and to conduct quantitative 
! decompositions.
!
! Refer to mod_sa.f90 for the simulated annealing 
! algorithm.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


    Program main

    use parameters
    use globals
    use tools

    implicit none
    include "mkl.fi"


    !CALL MPI_init(ierror)
    !CALL MPI_comm_rank(MPI_comm_world, proc_id, ierror)
    !CALL MPI_comm_size(MPI_comm_world, num_proc, ierror)

    !IF (proc_id==master) THEN
    !    t1=mpi_wtime()
    !    call cpu_time(t1)
    !end if



    if (comparastat==0)then
        call date_and_time(values=time1)
        diter0=time1(5)*3600+time1(6)*60 +time1(7)+time1(8)*0.001;
        call siman()

        call date_and_time(values=time1)
        diter1=time1(5)*3600+time1(6)*60 +time1(7)+time1(8)*0.001;

        print*,'Elapsed time--simulation',diter1-diter0
        ! CALL MPI_FINALIZE(iError)




    else
        call compstat()
        !CALL MPI_FINALIZE(iError)
    end if







    pause

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    contains

    subroutine siman()

    use parameters
    use globals
    use tools
    use simulated_anneal

    implicit none


    integer::nfcnev
    integer, parameter:: poe=poefirm, estimation=estimate
    logical::max
    real*8:: temper, rt,fopt
    real(8):: x0(nparam),x0read(nparam),lb(nparam),ub(nparam),vm(nparam),step(nparam),xopt(nparam),eps0 !,modelmom2(nmom,1),modelmom2_w(nmom,1)
    INTEGER     ::  nt, ier, iseed1, iseed2,  maxevl, iprint,  nacc, nobds,  neps

    integer :: date_time(8)
    character (len = 12) real_clock(3)
    character (len = 29) resfil
    character (len = 29) readfil
    character (len = 35) startfil
    character (len = 35) endfil
    character (len = 1) junkfill

    character (len=50) :: covfil
    !declarations
    integer :: numevals,ct,evalst,evalct,paramct,seeddim,numstatics,newrun,npart,psoseed,maxpsoit,xquicknum
    integer, allocatable :: seedarray(:)

    double precision :: thetavec(nparam),phipso(2),xpsotol,fpsotol,xquicktol

    
    newrun =2 
    psoseed=2501

    npart = n_proc

    xpsotol = 1.0e-3
    fpsotol = 1.0e-3
    xquicktol = 1.0e-3
    xquicknum = 2
    phipso = (/2.05,2.05/)
    maxpsoit = 5000


    max = .false.
    neps=4
    eps0 = 1.0D-5
    rt = .85
    iseed1 = 1
    iseed2 = 2

    nt = 5
    iprint = 2
    temper=5.0

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    maxevl =20000 
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!




    do i=1,nparam

        step(i)=0.275

    enddo



    vm=0.5

    if (manual == 1) then


        x0=(/ 0.64393  ,     3.6360   ,      0.69697,   0.51376E-01 , 58.194 ,     0.37432,      0.86264 ,     0.63294E-01   ,  0.24349E-01  /)


        lb=x0-0.3*x0
        ub=x0+0.3*x0





    else

        covfil = "Datfils/x0read2.txt"
        open (unit=62,file=trim(covfil),status="old")

        read(62,*) x0read

        close (unit=62)
        x0 = x0read

        !print *,'x0',x0

    endif




    CALL CPU_TIME(t1)




    xx=x0

    if (estimate==0)then
        PRINT *, "solving the model... "

        call fcn(nparam,x0,fopt,poe)
        OPEN (UNIT=1,FILE='params.txt',STATUS='replace')
        WRITE(1,FMT='(f20.15)') x0
        CLOSE (UNIT=1)

    else
        PRINT *, "Estimating the model using estimation method ", estimate_method
      
            call sa(poe,nparam,x0,max,rt,eps0,ns,nt,neps,maxevl,lb,ub,step,iprint,iseed1,iseed2,temper,vm,xopt,fopt,nacc,nfcnev,nobds,ier)
     
    end if








    CALL CPU_TIME(t2)


    PRINT *, "Total running time was ", t2-t1 ,"seconds"

    end subroutine siman



    subroutine compstat()

    use parameters
    use globals
    use tools

    implicit none


    real*8::allmoments(nmom,nparam,nvary),x0(nparam),fopt
    integer, parameter:: poe=poefirm
    double precision :: xstore(nparam,nvary)



    x0=(/ 0.64393  ,     3.6360   ,      0.69697,   0.51376E-01 , 58.194 ,     0.37432,      0.86264 ,     0.63294E-01   ,  0.24349E-01  /)


    comp_lb=x0-0.01*x0
    comp_ub=x0+0.01*x0


    maxmin(:,1)=comp_lb
    maxmin(:,2)=comp_ub

    OPEN (UNIT=1,FILE='minmax.txt',STATUS='replace')
    WRITE(1,FMT='(f20.15)') maxmin
    CLOSE (UNIT=1)




    do np=1,nparam !param loop
        do vary=1,nvary

            x0=(/ 0.64393  ,     3.6360   ,      0.69697,   0.51376E-01 , 58.194 ,     0.37432,      0.86264 ,     0.63294E-01   ,  0.24349E-01  /)

            !perturb parameters one by one, only change one parameter each time
            x0(np)=maxmin(np,1)+(maxmin(np,2)-maxmin(np,1))*(vary-1)/(nvary-1)


            xstore(:,vary)=x0(:)
            print *, 'parameter', np
            print *, 'vary', vary


            !record all moments for each parameter in a huge matrix
            call fcn(nparam,xstore(:,vary),fopt,poe)
            allmoments(:,np,vary)= modelmom_nvary(:)             !allmoments2((i-1)*nmom+1:(i-1)*nmom+nmom)
        end do
        


    end do

 

    OPEN (UNIT=1,FILE='allmoments.txt',STATUS='replace')
    WRITE(1,FMT='(f30.15)') allmoments
    CLOSE (UNIT=1)
 

    OPEN (UNIT=1,FILE='xstore.txt',STATUS='replace')

    WRITE(1,FMT='(f30.15)') xstore
    CLOSE (UNIT=1)



    end subroutine compstat


    end program main
